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ABSTRACT 

Characterization of transiting planets with transit timing variations (TTVs) requires understanding 
how to translate the observed TTVs into masses and orbital elements of the planets. This can be 
challenging in multi-planet transiting systems, but fortunately these systems tend to be nearly plane- 
parallel and low eccentricity. Here we present a novel derivation of analytic formulae for TTVs that are 
accurate to first order in the planet-star mass ratios and in the orbital eccentricities. These formulae 
are accurate in proximity to first order resonances, as well as away from resonance, and compare 
well with more computationally expensive N-body integrations in the low eccentricity, low mass-ratio 
regime when applied to simulated and to actual multi-transiting Kepler planet systems. We make 
code available for implementing these formulae. 

Subject headings: planets and satellites: detection 


1. INTRODUCTION 

No planet orbits on a precisely Keplerian orbit: post-Newtonian corrections, stellar oblateness, a nd, most impor¬ 
tantly, planetary perturbations cause deviations from a periodic ephemeris for transiting exoplanets (Miralda-Escude 


|2004[ [Holman fc Murra^|2005[ Agol et al. |2005| [Heyl fc Gladman|2007[ jlNesvorny fc Morbideili 
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Schneider! 
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2010). Transit-timing variations (T1 


blesvorny et ai. 


i 2al |Nie||2013t |Aie et a l. 

^ TN esvoriiy et al.||2^ 


7V s) have been used to confirm that transit signals are in tact 


due toTlanets (Holinan et al. 2010 Ford et al. 2012b Steffen et al. 2012 Steffen et al. 2013 Fabrycky et al. 


2014), to detect and characterize non-transitmg planets ([Ballard et al. 


2012 


mr 


JT5|7 and to r nake p recise measurements of the masses and dynamical states o: 


multi-transiting exoplanet systems (e.g. Carter et al. [2012 1. 

For the latter two applications, fast computation of 'IT vs is required for rapid searching through parameter space 
for perturbing companions, and for rapid computation of the posterior distributions of the masses and orbital elements 
of transiting planet systems. Numerical computation of TTVs can be sped up through symplectic i ntegration, through 
a more efficient numerical solution of Kepler’s equation, and through transit time interpolation (Deck et al. 2014); 
however, this approach can still be too computationally intensive for high multiplicity systems, and does not pinpoint 
the physical origin of constraints upon planetary system properties. Analytic formulae based on perturbation theory 
can greatly speed computation, but the perturbation theory to high order in eccentricity and inclination becomes 
comp licated quickly, and numerical codes that implement the analytic form ulae have not been released or widely 


used (Nesvorny & Morbideili 
order (in eccentricity) pertur 


2008 Nesvorny||2009 Nesvorny & Beauge 2010). Much can be accomplished with first 
aation theory because orbital eccentricities of many planets exhibiting TTVs are small. 


TTVs are most easily observ ed for pairs of planets near a mean motion resonance; thanks to the low eccentricity of 


the systems in consideration (Fabrycky et al. 2014 Hadden & Lithwick 
Albrecht||2015 1, the first order resonances are most represented among 


2014 Limbach & Turner||2014 Van Eylen & 

_ T'V pairs, and it is tor hrst order resonances 

that a hrst order theory is adequate. TTVs caused by these fir st order resonant interactions are primarily sinusoidal 
and are subject to a degeneracy between mass and eccentricity ( Bone et al.|2012 ), caused by mixing of two frequencies 
of perturbation which are aliased at the frequency of the transiting planet, as explained in an elegant analysis by 
Lithwick et al. (2012, hereafter LXW12). To break this degene racy requires the measurement of add itional modes, 
such as the short-timescale TTVs known as ‘chopping’ variations (iHolman et al.|2010' Deck & Agol 2015), or statistical 


analysis of many systems (Wu fc Lithwick|2013 


is only approximate, and breaks down tor pairs further from resonance (Deck fc Agol|2U15). These issues motivate the 


Hadden fc Lithwick|2014 ' Xie| 2014|. In ad dition, the LXW12 analysis 


current paper in which we derive an explicit formula for TTVs accurate to hrst order in eccentricity and planet-star 
mass ratio, valid for (nearly) plane-parallel transiting planets (although Nesvorny & Vokrouhlicky||2014 showed that 
mutual inclinations of planets can be large and still be well described by coplanar TTVs). 

We expect that these results will be useful a) f or determining how different frequencies within the TTV signal 
constrain the planetary masses and orbital elements (Deck fc Agol|2015|; b) for analyzing systems with a large number 


of interacting transiting planets by making linear additions ol the analytic formula for pairs of planets (Lissauer et al. 


agol@uw.edu 


















































































































2 


Agol & Deck 


2013 Jontof-Hutter et al.||2014|; and c) for rapid search through parameter space ofjaerturbing planets. 

We hrst summarize the 'i"i'V solution to first order in eccentricity and mass in ^ (the full derivation is given in 
appendix]^. We then compare these with prior results, both analytic and numeric ((|^, including comparison of 
analytic and numeric analyses of specific systems. We discuss the numerical implementation and speed in We end 
with a discussion of the possible applications and future directions (ij^. 

2. FIRST-ORDER SOLUTION: 

Here we give a complete summary of the assumptions and variables used, and the solution to the first-order equations 
for readers that wish to simply use the results of this computation. The details of the derivation are given in Appendix 
Since multi-planet transiting systems typically have nearly edge-on orbits and hence small mutual inclinations, it is 
usually sufficient in analytic approximations to treat the problem in the plane-parallel approximation. This leaves four 
orbital elements for each planet (semi-major axis, Oi, mean longitude, A^, eccentricity, e^, and longitude of periastron, 
Wi), plus the mass ratio of each planet to the star, rui/m*, where mi is the mass of the inner planet, m 2 is 

the mass of the outer planet, and m* is the mass of the star. For nearly circular planetary orbits, there are two small 
dimensionless parameters in the problem: = mi/m* and Ci. The usual procedure for computing transit timing 

variations is to: 1) to write down a Hamiltonian (or disturbing function) for perturbations due to another planet; 
2) expand the Hamiltonian as a function of the orbital elements to the order in eccentricity desired plus one (e.g. if 
a transit timing solution is needed to first order in eccentricity, then the Hamiltonian must be expanded to second 
order in eccentricity), including the linear combinations of mean-longitudes leading to the important resonant terms 
necessary for sufficient accuracy; 3) compute the variation in the orbital elements using Hamilton’s equations, which 
are four first-order partial differential equations for each planet, and involves differentiating the Hamiltonian with 
respect to the orbital elements (which can be a rather complex operation); 4) integrate the resulting equations as a 
function of time; 5) compute the true longitudes, 9i = as a function of time, where Oi^x is the unperturbed 

Keplerian orbit, and S9i is the perturbation of the ith planet caused by its planet companion(s); 6) compute the transit 
timing variations: 

6U = -9~x5Bi. ( 1 ) 


This is the approach taken by Agol et al. (2005), Nesvorny fc Vokrouhlicky (2014), and L XW12. A diffe r ent ap proach 
employing Hamiltonian perturbation theory (JNesvorny &: IVlorbidelli 2008) was used in Deck & Agol (2015). This 
involves determining the canonical transformation between the full canonical orbital element set and the average set; 
the TTVs, which are deviations from an average “Keplerian” orbit, can be derived from this transformation. 

The standard procedure outlined in detail above (based on Hamilton’s equations) has the advantages of requiring 
only first-order differential equations for the computation, and the advantage of using standard methods in celestial 
mechanics for the computation. However, there are two possible drawbacks: 1) the expansion of the Hamiltonian in 
orbital elements can be rather complex; 2) the main quantity of interest for transit-ti ming variations is 59j, rathe r 
than the perturbed orbital elements. The derivation based on canonical transformations ( Nesvorny fc Morbidelli|2008 ), 
though elegant, has the disadvantage of requiring the extra machinery and knowledge of Hamiltonian perturbation 
theory. 

In our new derivation we forgo computing the orbital elements, and simply treat the problem in polar coordinates 
(ri,9i). We then use Newton’s equations in terms of a disturbing function which can be expressed as a function 
of polar coordinates, with the added advantage that Newton’s equations make clearer which forces are causing the 
perturbations. This approach has some possible advantages: 1) only two differential equations are necessary (albeit 
second-order rather than first-order); 2) the derivatives of the disturbing function with respect to the polar coordinates 
are easy to compute; 3) the perturbed polar coordinates directly yield the transit timing variations; 4) the resulting 
expression is more compact than in the Hamiltonian formulation. The second-order differential equation may seem 
like a drawback, but it can be solved using complex notation (as in LXW12) and by expanding the derivatives of the 
disturbing function in terms of orbital elements, which yields harmonic functions which are easy to integrate. T he final 
answer is expressed as a sum over harmonics of the perturbing planet’s orbital frequency (Deck & Agol 2015). Each 
coefficient for each planet in the harmonic series solution can be solved for by inverting three two-by-two matrices, 
which have a standard format, resulting directly in the transit timing variations at a particular frequency. 

The unperturbed orbital frequencies, Ui = 2TTjPi, are defined by n? = Gmi,/af. As usual, a = 01/02 « {Pi/■ 


We define Aj^n = a 


n+l a" 


1 2 daY^da2 


(02 


where is the Laplace coefficient. 
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r2n 


dO- 


cos {j9) 


{1 + — 2a cos 0)1/2 

The difference in mean longitude of the planets is '0 = Ai — A 2 . Auxiliary dimensionless quantities are: 

/3j = j(«i - n2)/ni = j{l - a^/^), 

Kj = j{ni - n 2 )/n 2 = j{a~^^'^ - !)■ 

The functions Ajmn we use below are given by: 

AjOO = ^ 1/2 


( 2 ) 


(3) 
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Table 1 

Coefficients for tx( 7 ,ci,C 2 ) and v±{(^,di,d 2 ) in first-order TTV solution (the v± coefficients correspond to the k values in brackets [..]). 
i ±k 7[C]ci[di] C2 [d2\ 


1 

0[±1] 


“i ( 

j 

a (Ajio - OiSji^ 

1 

±1 


“j ( 

ijAjoo — |Ajio + |(i 7 2)o(5ji^ 

a (ijAjio — §Aj20 7 

1 

±2 

dj ± 

“i ( 

^j-^jOO ~ ~ (1 T 

Of (^j^jlO 2 -^jll (1 ^ 

2 

0[±2] 

Kj 

-3 1 

f Ajoo — 

Ajoi — 

2 

±1 

K,j di 

-3 1 

(ijAjoo — lAjio — (1 ± l)a~'^Sjij 

±iAjoi — |Ajii — (1 ± l)a~'^Sji 

2 

±2 

Kj ± 1 

-3 1 

(TiAjoo — + 5(1 ± 2)«“^<5ji^ 

7iAj01 — ^Aj 02 ± 


AjiQ = adb^^j^/da, 

Aj2o = c? 9 ^ 61/2 / ’ 

AjQi = —{AjiQ + AjQo) = —{adb^^j^/da + 

Aj 02 = 2AjQQ + AAjiQ + Aj 2 o = ‘ 2‘^^2 da + a^d'^b^^^^jdo?^ 

Ajw =-{2Ajio + Aj2o) = -2adb^fj^lda - c?d'^b^[j^ldc?. (4) 

To use this solution in computing TTVs, the longitudes need to be computed from the observed transit times; the 
mean ephemeris, (to^i,Pi), may be used in computing the (unperturbed) orbital ephemeris. Now, the mean longitudes 
are given to first order in eccentricity by 


Ai = 27r 


/ i — to,i 


+ 2ei sin-nui. 


(5) 


if we assume that the orbital reference is along the line of sight, and thus Ai « 0 at the times of transit. 
The solutions for the inner planet (z = 1) and outer planet {i = 2) are given by: 


[fiJ Ibii’ - (Ai - wi)] + /}sin [jz/; + (Ai - zui)] 

^ i>i 

+ /iy-\e2 sin [jii - (Ai - -njs)] + /|j+\e2 sin [jip + (Ai - 1x72)] , 

St2 = [■^ 2 °- + ^2 sin [j^|; - (Aa - wa)] + /i^^^ea sin [j^ + (Aa - wa)] 

i>i 

+ /Ij+iei sin [jtp - (Aa - wi)] + /a^i^ei sin [jz/> + (Aa - wi)] , 


( 6 ) 


where the functions are given by: 


/ij (“) = '“( 7 , Cl, ca) + S^kV±{C, di,d2), 
(3 + 7^) Cl + 27C2 


'* j 

m(7,Ci,C2) = 


C±(C,rfl,C^ 2 ) = - 


(7) 


( 1 -^ 2 ) 

(± (1 — C^) + 6C) di + (2 + (C^) da 

C(i-0(C±i)(C±2) 

where 7 , ci, and ca and C, di, and da are given in Tabled and Sik is the Kronecker delta function. Note that the top 
signs in ±,=F correspond to +k values, while the bottom correspond to —k. The functions fP are solely a function 
of j, ±fc, and a. 

In practice the sum over j from 1 to 00 must be truncated at a finite value of i mnx • Typically im.ax does not need to 
be chosen to be too large since the Laplace coefficients decline in amplitude with j (Deck & Agol 2015). We recommend 
choosing a j, 


large enough such that the resulting computation is converged. 

As an example of using Table ll the coefficient has z = fc = 1, 7 = /3j — 1, Ci = aj (— 


and C 2 = a (—jA 
is given by: 


(—jAjoo 


~ + ^aSji 


jio - 2 ^i 20 


+ aSji^, = Pj, di = aj (Ajoo — adji j, and d 2 = a (a^iq — aSji^. Then, the coefficient 


/v 


( ^^(a) = u{/3j - l,ci,C 2 ) + z;_(^j,di,d 2 ). 


( 8 ) 


3. COMPARISON WITH OTHER FORMULAE 
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Figure 1. Dimensionless coefficients, ^ for the inner planet. The dashed lines show where the coefficients are negative. 


The zeroth-order solution (in t he lim it e^ = 62 = 0 ) comp ares exactly with the results given in Agol et al. (2005), 
Nesvorny & Vokrouhlicky (2014), and Deck & Agol (2015) which were derived with Hamilton’s equations and the 
approach based on canonical transformations; this is reassuring given the very different approach used in this derivation. 
We have als o rederived the first- order eccentricity equations using the approach based on canonical perturbation theory 
employed in De ck fc A^ol (2015), and found exact agreement with the results presented here to first order in eccentricity. 

In Figures^^ and (j2|we plot the eccentricity-dependent co efficients, as a function of period ratio, P 2 /P 1 « 


v-3/2 


a . The zeroth-order eccentricity coefficients are plotted in Deck & Agol (20151. The first-order coefficients can show 


three singularities for the terms with superscripts (— 1 ) and (— 2 ) near hrst-order resonance, second-order resonance, 
and a = 1 . 


3.1. Comparison with first order resonant equations 

LXW12 present a formula valid near (but not in) first order mean-motion resonances that captures the behavior of 
resonant terms in an elegant, but approximate, manner. Here we compare the complete formulae given here to their 
near-resonant formulae. 

The expressions for u and v± do not show the same dependence in the denominator as the expressions in LXW12; 
their expression just contains the resonant frequency, jni — {j + l)n 2 , while ours contains additional frequencies. We 
carried out the partial fraction expansion of u to isolate the denominator which matches LXW12’s expression, and we 
find that the expressions agree exactly with their expressions (A28) and (A29). 

To compare our full expression with LXW12’s, we have computed the eccentricity-dependent j = 2 (near 2:3) 
expression for the inner and outer planets (as this term is unaffected by the indirect terms). We re-write the TTV 
formulae derived here and in LXW12 for the f—th planet perturbed by planet k as 


5ti = 




E 


APe, sin (jA* -h (j)P) + Alj^eu sin (jAfe -h 


(9) 


where we have set = Ai = 0 at transit (this incurs some error, at order e, but that is a second order effect since we 
are comparing the TTV term linear in e). When written in this way, the amplitude and phase depend only on a^vji, 
and 1372. 

For the 3:2 resonance, the LXW12 resonant term depends on ei/(2ni — 3 ^ 2 )^ and e2/(2ni — 3 ?t. 2)^, which both decline 
quickly away from resonance, and thus other terms that depend on ei and 62 make a more significant contribution 
further from resonance; hence our formulae agree close to resonance but diverge away from exact commensurability. 
Figurel^shows the fractional error in the amplitude A and phase (j) of the terms that are proportional to the eccentricities 
of the Janets for w\ « -072 = 0.45 radians. The error in the LXW12 expression (A28) and (A29) reaches ^20% at 
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Outer planet 



fife) 

Figure 2. Dimensionless coefficients, ^ ' for the outer planet. The dashed lines show where the coefficients are negative. 

a 5% separation from exact resonance in this case; if we use the further approximate expression given in their main 
paper in lieu of (A28), the discrepancy for the inner planet in creases to ^30%. This is similar to the error in the 
zeroth-order component of their expression ( Deck fc Agol||[^015 ). The zeroth- and first-order eccentricity terms have 
a different dependence on longitudes: for example, for the inner planet the zeroth order term scales as 
while the first order term scales as where i = v^—T. Since the mean longitude of the inner planet is 

nearly identical at each transit of the inner planet, the Ai term in the exponent is approximately constant, while the 
A 2 dependence is identical in both terms, leading to aliasing of these coefficients. Thus the error incurred in their 
approximation can lead to a different phase dependence and amplitude for this aliased term away from resonance. 

3.2. Comparison with N-body integrations 

We have carried out extensive integrations of three-body_systems using TTVFast (|Deck et al.||2014|), and compared 
the results with the first-order analytic formulae (equation]^. Note that the TTVFast code uses the convention of the 
longitude of periastron being measured from the sky plane to match the convention of ra dial velocity surveys, while 


here we use the observer’s line of sight as the reference direction, as done in LXW12 and [Deck & Agol| (2015). The 
longitudes computed from the TTVFast code need to have 7 r /2 subtracted to make the plots shown below. In addition, 
the orbital elements accepted by TTVFast are the instantaneous/osculating orbital elements (initial conditions) at the 
specified initial time, while the orbital elements used in these formulae are the mean orbital elements of the planets 
over the timescale of the observations. 


3.2.1. Eccentricity and period ratio dependence 

Figure ^ compares the precision of the analytic formula as a function of a = {Pi/P 2 )‘^^^ and eccentricity of both 
planets, ■^ich are set to be equal, ei = 62 . We have set vui = + tt, which we found (approximately) maximizes the 

discrepancy of the analytic model compared with the N-body model, and vj\ = vji which (approximately) minimizes 
the discrepancy; hence the figures bracket the precision of the analytic model. This is due to the fact that the anti¬ 
aligned longitude geometry causes the planets to be closer at conjunctions that occur when the inner planet is at 
apoapse and the outer is at periapse; their proximity at these conjunctions causes their gravitational interactions to be 
more sensitive to deviations from the epicyclic approximation, which are second order in eccentricity, and thus missing 
from our computation. We have assumed that the period of the inner planet is Pi = 30 days, and we have integrated 
the system with TTVFast for 1600 days, about the duration of the initial Kepler mission, assuming plane-parallel 
orbits. For these tests we assume tn = h -2 = 10“®, and we selected random values for the longitudes of the planets 
at the initial time. For each set of initial conditions, we output the orbital elements at regular intervals during the 
N-body integration, from which we computed the average orbital elements over the duration of the integration. These 
averaged orbital elements were used for computing the amplitudes of the analytic model, which we summed up to 
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Figure 3. Fractional error in the coefficients of the 3:2 resonant TTV expression given in LXW12 compared with the j = 2 terms of our 
analytic expression. Black indicates dependence on the ei term, while red indicates the dependence on the 62 term. The top show the 
fractional errors in the amplitudes of the inner (left) and outer (right) planets. The bottom shows the fractional error on the phases. 

jmax = 10. We optimized the fit of the analytic formula to the numerical TTVs by allowing the ephemerides of the 
planets to vary in the formula, but holding the eccentricity vectors and mass ratios fixed at the values computed from 
the time-averaged N-body simulation, while we computed a in the analytic formula from the ratio of the periods 
derived from the best-fit ephemerides. 

The fractional precision was computed from the RMS of the residuals of the best analytic model ht to the TTVs, 
divided by the RMS of the TTVs computed from the N-body integration. Figure]^ shows that the formula works to 
better than 10% precision for a wide range of a—e parameter space. However, it fails near resonances, most signihcantly 
for the j'.j + 1 resonances indicated in green, and j:j -|- 2 in blue. For the outer planet, there are narrow regions near 
l:j period ratios for which the formula does poorly. The disagreement grows in breadth for larger eccentricities. This 
diagram can be used to pinpoint the relevance of the analytic formula for a particular system, and we suggest that the 
analytic formulae should be used with caution in the regions where the formula disagrees by more than 10% precision. 

Most of the regions where the formula fails are near resonance. In these cases, the residuals can frequently be fit 
by sinusoidal variations at the relevant resonant frequencies of the higher order resonant terms that are not captured 
in the first-order model; when including these sinusoidal terms in the fit, the residuals drop dramatically near the 
resonances. Thus, the analytic first-order solution plus a sinusoid with arbitrary amplitude and phase can be used for 
systems in which only the shape of the transit timing variations plus the specific variations of the non-resonant terms 
is necessary (although this approach breaks down for large enough eccentricity). 

3.2.2. Mass dependence 

We have carried out simulations for a range of masses, keeping mi = m 2 (/ii = pi^)- We find that the fractional error 
of the analytic formula grows near resonances and near a = 1 as the mass increases with weaker dependence upon 
eccentricity. For small a, the formula works well up to mi = m 2 = 10“^m*, while near 1:2 period ratio, for example, 
the error broadens around the resonance. 

Figurej^shows the fractional error in the formula (computed as in the eccentricity dependent case) for Ci = 62 = 0.001 
(the mass dependence of the precision is nearly independent of eccentricity) with w\ = -072 -I- tt (the results look very 
similar for vji = W 2 )- The formula is accurate for a broad range of masses, but for some systems, such as Planet 
Hunters 3c/d, indicated in Cyan in Fig. the discrepancy becomes large, « 10% (the masses and eccentricity vectors 
for this plot differ from PH3c/d, but a juot made for the parameters of that pair of planets looks very similar). 
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Figure 4. Fractional precision of the analytic formula compared with TTVFast. Left: Aligned longitudes of periastron {vji = vj 2 )'i Right: 
Anti-aligned longitudes of periastron (vji = + tt). The dotted lines indicate the 10% precision level. Cyan dots show the approximate 

position of Kepler-18b/c at the 85.15% posterior eccentricity value, while magenta is used for Kepler-28. The region in the upper right is 
Hill unstable; these models were not computed, and default to 100% uncertainty in this plot. The green dashed lines show the locations of 
j\j -h 1 resonances, while the blue dashed lines show y.j -h 2. 



0.2 0.4 0.6 0.8 

alpha 


Figure 5. Fractional precision of the analytic formula compared with TTVFast versus a and mass ratio of the planets to the star 
(mi = m 2 ). Top: inner planet; bottom: outer planet. The dotted line indicates the 10% precision level. Cyan dots show the approximate 
position of PH3 (although note that PH3 does not have equal masses; however, the plot is similar for parameters appropriate for PH3). 
The region in the upper right is Hill unstable; these models were not computed, and default to 100% uncertainty in this plot. The green 
dashed lines show the locations of j\j 1 resonances, while the blue dashed lines show j:j -h 2. 


3.3. Comparison with two-planet systems 

We have carried out fits to systems with two interacting planets described in LXW12, and we have re-fit Planet 
Hunters (PH3) c/d. We have carried out N-body dynamical analyses using TTVFast in addition to fits with the 
analytic first-order formula in order to assess its utility in analyzing multi-pla net systems. 

Our first case is Kepler-18c/d, which was publis hed by Cochran et al. (2011) and also analyzed by LXW12. We used 
the same transit times and uncertainties from the Cochran et al. 12011) paper to allow for direct comparison to their 
results; these transit times were also used in LXW12. We carried out a markov chain monte carlo (MCMC) analysis 
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Mass e cos(omega) 

Figure 6. Comparison of numerical and analytic analysis of the transit times of Kepler-18. Left: histogram of the masses of each planet. 
Right: comparison of the 68% confidence distribution of the eccentricity vectors of the two planets (18c solid; 18d dashed). In each case 
black indicates the results from the N-body analysis, while blue indicates the results of the analytic formula. 


2010 ). We allow a multiplic ative factor for the 
le eccentricity of each planet (Ford|2006). Figure 


using an affine-invariant population approach (Goodman & Weare 
timing uncertainties on each planet and we place a uniform prior on t 
shows a comparison of the results of the MCMC analyses with N-body integration versus the analytic formulae with 
jmax = 5. For this system the mass ratios are « 5 x 10“® and the eccentricities are of order 0 — 0.02 (1-cr), while 
a « 0.64, which is a regime in which the first-order formula is accurate to <9% compared with N-body integration (see 
Figure]^ the cyan dot indicates the approximate location of the upper end of the eccentricities of Kepler-18c/d). For 
conversion of the mass ratios to planet masses, we assume that m* = O.972M0 (we ignore the uncertainty on this stellar 
mass). The masses of the planets derived from N-body are mi{nbody) = and m 2 {nbody) = 

while from the analytic formula are mi(analytic) = 13.2'^^' qM^ and m 2 {analytic) = These are well 

within 1 — (T of one another, the mar kov chain posteriors s how very similar distributions (Fig. m. We note that these 
results differ from those reported by Cochran et al. (2011) who used N-body to estimate the transit times, found the 
best fit using Levenberg-Marquardt optimization, and estimated the uncertainties from the Hessian matrix at the best 
fit parameters rather than a full posterior analysis. Our results also differ from the estimates in LXW12, which only 
solve for a ‘nominal’ mass assuming zero eccentricity using the approximate near-resonant formula. We feel that our 
results should be superior to these prior results, and warrant a more extensive analysis with the full Kepler dataset as 
well as radial velocity measurements. 


We next compared analyses of the Kepler-28 system, which was originall y studied by S teffen et al. (2012), and 
included in the analysis of LXW12. AW used the transit times published by Steffen et al. 1 2012), and followed the 
same procedure as Kepler-18. Figure shows a comparison of the masses and eccentricity vectors for this system, 
which has a mean period ratio of {P 2 /P 1 ) = 1.52, just wide of 2:3 period ratio, corresponding to a = 0.7563. 
The masses are poorly constrained due to the degeneracy with eccentricity, which allows the eccentricity to wander 
to larger values. For conversion of the mass ratios to planet masses, we assume that m* = 0.89Mq (we ignore 
the uncertainty on this stellar mass). A comparison between the mass constraints from N-body and the analytic 
formula gives: mi{nbody) = versus mi{analytic) = and m 2 (nbody) = 5 . 11 ' 13 qM 0 versus 

m 2 {analytic) = 4 . 1 + 2 ^ 3 ^ 70 . The 85% confidence value of ei is 0.098, while for 62 is 0.076, while the longitudes of 
periastron within the posterior distribution are primarily anti-aligned; the location of these points is indicated with a 
magenta datapoint in Figure® Note that in the anti-aligned w case the analytic formula is valid to larger eccentricities, 
and thus is adequate to describe this system. 

As F igure ® indicates, the Planet Hunters 3 (PH3) system, the outer two planets (c/d) analyzed in|Deck & Agol 
(|2015|), has amr ge discrepancy due t o the large mass of the outer planet. The excellent agreement with the chopping 
formula given in Deck fc Agol] (2015|) is still imperfect; the figure in that paper was mistakenly produced with larger 


Schmitt et al. 


timing error bars than used in 
first-order for mula indicates. Figure 8 shows tbe results of a comparison of N-body and analytic fits to the PH3 transit 


(2014) which caused the agreement to appear slightly better than the 


times given in Schmitt et al. (2014); the planet masses assume to* = IMq. The analytic formula gives a significant 
discrepancy due to the large mass of the outer planet and due to the proximity to the 1:2 period ratio. 

To confirm that the large mass of PH3d led to this discrepancy, we took the best-fit parameters resulting from our 
N-body analysis of the real PH3 data, and reduced the masses of the outer two planets by a factor of 10. Using the 
outer two planets alone, we simulated transit times and added Gaussian noise at a level of 1/10 that of the noise of the 
real data to maintain the same signal-to-noise ratio as the actual data. We then modeled this simulated data using 
TTVFast and the analytic formulae. We found that the agreement between the N-body and analytic analyses becomes 
excellent (Fig. |^. 


3.4. Comparison with multi-planet systems 
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Figure 7. Comparison of numerical and analytic analysis of the transit times of Kepler-28. Left: histogram of the masses of each planet. 
Right: comparison of the 68% confidence distribution of the eccentricity vectors of the two planets (28b solid; 28c dashed). In each case 
black indicates the results from the N-body analysis, while blue indicates the results of the analytic formula. 




0.0 0.2 0.4 0.6 0.8 1.0 
M. 


Figure 8. Comparison of numerical and analytic analysis of the transit times of Planet Hunters 3c/d. Left: 68% confidence contour of 
the masses of both planets. Right: comparison of the 68% confidence limits with the mass of the both planets reduced by a factor of 10. 

The two-bod y solution can be use d for more than two bodies by addition of two-body TTV solutions for each pair 
of two planets (Lissauer et ah 11 ^ : 

( 10 ) 


^ 




where are the solutions from equation for the iiih planet due to the 12 th planet. The sum over j for each 

pair of planets can be carried up to jmax to give sufficient precision for that pair of planets that is smaller than the 
measured timing precision. 

We 


used the transit times reported in 
analytic TTV signal given as the 


Our first system of study is Ke pler-51 (Masud a 2014), consisting of planets with period ratios close to 1:2:3. 

Masudall 2014 pto carry out dynamical models with N-body /TTVFast and with an 

We included up to j = 6 


sum ot the IT'Vs of the three adjacent pairs of planets, 
in the TTV signals. The results show excellent agreement; Figure shows the measured transit-timing variations, as 
well as the best-fit N-body and analytic TTVs. The results of the MCMC analyses are compared in Figure [TO] which 
shows the posterior distribution of masses and eccentricities measured with both analyses, assuming m* = IMq. For 
two of the planets the eccentricities are consistent with zero; in these cases the tail of the eccentricity histogram is 
heavier for the N-body than the for the analytic formula. The masses of the inner two p lanets from our N-body and 
analytic MCMC analyses agree well with the masses from the analysis in Masuda (2014), while the mass of the outer 
planet is small by « 1 — cr. 


4. NUMERICAL IMPLEMENTATION AND SPEED 

The primary computational burden of equation (H lies in computing the Laplace coefficients. We use a series 
solution for these coefficients, which gives both speed and accuracy, using code shared by Jack Wisdom. The secondary 
computational burden is in computing the sine and cosine terms, which involves four angles, and thus requires eight 
evaluations. We carry out the computation of higher j sines and cosines using trigonometric addition formulae, which 
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Figure 9. Transit timing variations from |Masuda| if 
N-body model computed with TTVFast ( black J, ancf 
(blue). 


2014| 

the 


I for Kepler-51 (red; 0, 1, 2 stand for 51b,c, 620.02), compared with the best-fit 
oest-fit two-planet, first-order eccentricity formula summed over pairs of planets 


means we only need to compute eight trigonometric functions at each transit time once; the rest are gotten from 
addition and multiplication of these. 

In the cases that we run a Markov chain for a set of planets, the initial value of a is known fairly well from the period 
ratio of the planets. In this case the Laplace coefficients and their derivatives needed for the solution can be Taylor 
expanded at the a given by an initial fit to the transit times, and these coefficients can be stored for evaluation of the 
coefficients at slightly different values of a encountered during the MCMC simulation. This approach would not work 
if a is being varied over a grid (for example, in the case of searching for a perturbing planet with unknown period); 
however, computational effici ency can still be achieved by reusing the Laplace coefficients at different eccentricities 
(Nesvorny & Morbidelli||2008 ). _ _ 

We have coded the hrst-order formula, equation (|^, in C, IDL, Python, and Julia ( Bezanzon et al.|2012 ). We carried 
out a benchmark comparison of the Julia implementation of the formula with the C implementation of iTVFast, and 
we find that it is 400x faster when the Laplace coefficients are approximated from a Taylor expansion. As TTVFast 
is about 20 times faster than TTVs computed w ith standard N-body integrato rs, this represents nearly four orders of 
magnitude in speed up, similar to that found by Nesvorny & Morbidelli (20081. Note that if integer period ratios are 
chosen, sometimes the denominators of u and v± can become inhnite, causing divergence; we expect that this will not 
be encountered in practice as the formulae only apply to non-resonant planets. 

The code implementing these equations may be accessed at https://github.com/ericagol/TTVFaster 

5. DISCUSSION AND CONCLUSIONS 

In modeling transit timing variations, degeneracies and computational speed can each prohibit the accurate mea¬ 
surement of transiting planet masses and orbital properties, with their attendant uncertainties. The degeneracy due to 
aliasing near first order resonances (LXW12) can be broken with very high signal-to-noise due to the slight difference 
in the eccentricity d ependence as a function of period ratio, as well the presence of perturbations at other frequencies 
(Deck & Agol 2015). Here we have tried to improve the modeling of the terms which are linearly dependent upon 


eccentricity to provide a higher-fidelity analytic model to address both the degeneracy and the computation barriers. 

To this end, we have presented a first-order solution in eccentricity and mass ratio to the plane-parallel, near¬ 
circular 3-body problem on timescales shorter than the secular timescale. This improves to first order in eccentricity 
the original solution given in | Agol et al. (20051 which was der ived to zeroth order in eccentrici ty a nd hrst order in mas s 
ratio (this solution has also been given in different forms in Nesvorny & Vokrouhlicky 2014 and Deck & Agol||2015). 
The expressions are accurate compared with numerical integration over a wide range of parameter space relevant to 
the hundreds of multi-transiting planetary systems being found at short orbital periods with Kepler (Rowe et al.|2014 
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Figure 10. Comparison of numerical and analytic analyses of the transit times of Kepler-51 (0, 1, 2 stand for 51b,c, 620.02). Left: 
histograms of the masses of each planet. Right: histograms of the eccentricities of the planets. 


Lissauer et al. 
given in LXW 


20141. We find that this expression is more accurate than the stripped-down near-resonant formula 
2 , although their formula has a simpler form which clearly highlights the mass-eccentricity degeneracy. 


The first-order eccentricity formulae also can be used to model more than two planets with linear combinations of 
2-planet formulae, and works well for the system we tested here, Kepler-51. 

We used an approach starting with the Newtonian equations of motion rather than the Hamiltonian, and compute 
the perturbed polar coordinates of the planets’ orbits; this approach is akin to solving dispersion relations of differential 


Chandrasekhar 

1961 

Balbus & Hawley 

1998 


The unperturbed solution to these differential equations represents 


inhomogeneous component to the solution, which cause TTVs to vary at frequencies which depend upon integer 
combinations of the orbital frequencies of the two planets. Since the answer obtained in the end is the same as in 
the Hamiltonian (for C>(e°)) and canonical transformation approaches, this approach might be useful pedagogically 
for those more familiar with stability analyses. In addition, this approach might be useful for other problems, such as 
a stability analysis of a two-planet system or for carrying out the TTV computation to second-order in mass ratios 
Hi = mi/m*. The latter is interesting as it would reveal how TTVs of a transiting planet may be used to measure the 
mass of that planet (and not just the mass of the perturbing planet). 

We expect that t hese formulae will be used in carrying out initial fits to multi-transiting planets which show TTVs 
(Mazeh et al. J2^, in searching for companion perturbing planets to isolated planets showing TTVs, in characterizing 
multi-planet systems with TTVs to confirm and check for convergence of N-body MCMC analyses, in forecasting 
TTV amplitude for follow-up measurement, in estimating the optimum times for transit observation, and in making 
predictions for transit times to plan observations. It should be useful for estimating the densities, masses, and radii of 
the host stars and their exoplan ets (Agol et al.|2005 Montet & Johnson|2012 Hadden & Lithwick|2014 Kipping et al. 
'2014 Jontof-Hutter et al.||2015 ), and in comparing the IT'V solutions to radial velocity solutions tor t le masses and 


orbits of exoplanets. 'I'he analytic nature of our solution should be amenable to automatic differentiation (Fournier et al. 


20121, which could speed up optimizatio n based on gradient computation, and could also enable Hamiltonian/Hybrid 


Markov Chain Monte Carlo (Neal||2011 ). When a large number of planets transit a star and each show evidence for 
dynamical interactions, the number of free parameters describing the system becomes large, and thus MCMC becomes 
prohibitively computationally expensive. The first-order analytic formula developed here can be used for modeling 
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these systems if their masses and eccentricities are in the allowable range. As the formula is about 400 times faster 
than TTVFast, which is already about 20 times faster than Bulirsch-Stoer based integrators, the total speedup of 
about 8000 should make running chains long enough to conver ge more feasible, especially in tandem with parallel 


computation which can be easily adapted for population MCMC (Foreman-Mackey et al. 


2013). The analytic formula 


also has the advantage of being able to pinpoint which features ih the IT'Vs con strain the paramete rs of the system 
(LWX12). Since the TTVs of a planet display harmonics of the perturbing planet (Deck & Agol|2015), the amplitudes 
and phases of each of these harmonics can be measured directly from the TTVs, and then these can be used to place 
individual constraints upon the masses and eccentricity vectors of the planets. The regions where the co nstraints 


overlap ma y reveal the consistency and uniqueness of the solution for the system parameters in some cases (Deck & 


Agol 20151. 


Although in principle TTVs allow for unique measurements of planet mass and eccentricities, degeneracies between 
these parameters are often found for systems with low signal to noise. In these cases, the eccentricities can become 
extremely large, indicating unstable orbits, as long as the masses are adjusted in a corresponding manner. We applied 
Hill stability to avoid this problem when using the analytic formulae; the full N-body computation avoids this issue 
naturally since large eccentricities introduce second-order (and higher) variations (that our calculation ignores) which 
prevent the high-ecc entricity cases from f i tting the data well. It may be possible to break some degeneracies with transit 


duration variations (Pal & Kocsis 2008 Nesvorny et al. 2013), which can be computed with the same formalism we 


have described here, albeit in the plane-parallel limit. 

The first-order formulae described here could be extended to higher order in eccentricity and/or mass ratio, albeit 
with much more computational effort. A slightly more accurate formula might be obtained by computing the longitudes 
from Kepler’s equation at the times of transit of each planet rather than using the first-order eccentricity formula ([^, 


as well as using the exact formula for 9i at the transit times in equation ([^; this requires very little additional 
computational effort, but SOi will still be only accurate to first order in eccentricity. The solutions for the perturbed 
polar coordinates, {6ri,69i), can be derived in the same manner that we have derived the transit timing variations. 
These are needed for carrying out the higher-order perturbation solutions, and in turn could be used for modeling 
astrometric variations, radial velocity vacations, and pulsar timing variations of host stars to account for the interactions 
of planets to first order in eccentricity. 
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APPENDIX A 


A NEW APPROACH TO TTVS 

Here we give the detailed derivation of the first-order solution presented in Section 

5.1. TTVs from angular and radial variations 

We start with the equations of motion in Murray & Dermott (MD), 6.10-6.11, which are in heliocentric coordinates. 
We use the index i to label the planets, and denote the inner planet with i = 1 and outer with i = 2. TTVs can 
be computed from the true longitudes, 9i and 02 , with respect to the star; hence heliocentric coordinates are ideal 
for transit timing computation. We convert the equations of motion to polar coordinates, and isolate the radial and 
longitudinal equations: 


d 


on. 


+ “^rifiOi — -?^{Ui + TZi) — — li, 


09 ,. 

0 


'I’i ~ fi9j^ — (Ui -\- TZi) ■ 


( 11 ) 


The first equation can be rewritten as the time derivative of specific angular momentum, li, and without the disturbing 
force it expresses conservation of angular momentum, while the second equation includes centripetal acceleration in 
the radial direction as the second term on the left hand side. The term Ui = G{mi, mi)lri is the standard 
Keplerian potential, while TZi is the disturbing function which reflects the gravitational potential energy of planet- 
planet interactions. Because the plane t-pl anet interactions lead to only small perturbations of the base Keplerian 
orbits, we seek a solution to equations (11) of the form: = ri ^ -\- Sri and 9i = 9iji -|- S9i, where 9i j() is the 

unperturbed Keplerian polar coordinates of the orbits, and {Sri,69i) are the small perturbations. We can then plug 
these solutions into the equations of motion ( |11[ ) and expand in powers of Svi, 69i, mass ratio, and eccentricity. 

We normalize by Oi, which is the semi-major axis of the unperturbed Keplerian orbit (we do not perturb Oi, zui, or 
Bi, so they are fixed at the mean, unperturbed values). Then ri/at = l+Ci, where Ci is a dimensionless radial coordinate, 
so that ri/oi = ki and Srijoi = Sti- The solution to the unperturbed orbits to first order in eccentricity is 9i^K = 
Xi — Re [2iiz*e’'^’] = Ai -|- 2ei sin (A^ — vJi) and = —Re [z*e*^’] where i = VM, Zi = = ei{coswi -|- isinnji) 

is the complex eccentricity vector (LXW12), 2 * is the complex conjugate of z, and Re[.] is the real part. 


5.2. Perturbed equations of motion 


The perturbed equations become: 


...... . QTZ 

rjKS^i = -‘^fi,K9i,KSri - 2ri^K9i,KSri - 2ri^K9i,KSri - 2ri^Kri,K59i + 

Uwi 

r- o 5 i:5 , 52 r , 2 G(to*- b m*) ^ , &Ri 

Sri = 2ri^K9i,KS9i -b 9^ j^Sr,, + 3 Sri + ——. 


(12) 
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In these equations, we have cancelled the terms for the unperturbed Keplerian on both sides of the equation. From 
hereon we will drop the ‘AT’ subscript from the unperturbed Keplerian orbital elements. To first order in eccentricity, 
the TTVs of the tth planet are: 


SU = -n-\Se, + 2eM°^), 


(13) 


where is the perturbed solution to zeroth order in eccentricity. To lowest order in eccentricity, e* 
so: 


€-i COS (Aj ZUj), 


(14) 


5 . 3 . First order in eccentricity equations 

Let the specific angular momentum equal U = rfOi. Since U is a constant, 9i = —2lifi/rf. To first order in 
eccentricity, U = Uiof, where rzi = 2TrjPi. Substituting these into the above equations, and expanding to first order in 
eccentricity, we find 


dh - 2ni{l - ei)59i - - I0ei)6ei = 

(1 + 2,ei)59i — 2niei6ei + 2ni{l — ei)Sei + 2ei59i = qq. ' 


(15) 


Note that the terms with ei or q are first order in eccentricity; hence, the other quantities in these terms need to 
only be expanded to zeroth order in eccentricity. Denoting the zeroth-order solutions as and 59'f\ we find the 
differential equations governing the transit timing solution: 


(' 


69i + 2ni5ei = -2 ( eid'df^ + + 2ni {iiSef"’ 

1 dn^ 


lide. 


(0) 


1 dn. 


dii — Sn^Sci — 2ni69i = —2niei (^69^^^ + + 


(16) 


We assume that both sides of these equations are complex, and the final solution is found from taking their real parts. 

The quantity TZi is the disturbing function, which for the inner planet can be broken into two pieces: TZi = 
^^^{TZd + (xTZe) (MD 6.44), where TZ-d = a 2 /|r 2 — ri| and TZe = — (ci/ai)(a 2 /?’ 2 )^ cos {9i — 02 )- For the outer planet, 
there are also two pieces, 7^2 = + 0 (~^TZi), where TZj = — (1 + e 2 )(l + ei)“^ cos {9i — 6 * 2 ) (MD 6.45). As usual, 

a = 01/02 « {Pi/P2Y^^- 


5 . 4 . Expansion of the disturbing funetions 

The expansion for TZd is given in MD 6 . 66 . We are considering the plane-parallel case, so 4* = cosif — cos {9i — ^ 2 ) = 
0, in which case we only need to include the oc 4**^ term. Also, we would like a solution that is first order in eccentricity 
(of the unperturbed Keplerian orbit), so the term in brackets in MD 6.66 needs to be expanded to second order in 
Ci = — 1 (noting again that Ci is first order in eccentricity): 


i 


1^0 ’ k^O 


e'l'ei) "Aoj^k,i-k — Aoj-, 0,0 + Aoj-p^oci + Aqj, 0 , 1^2 + 2 ^ 00 . 0,262 + Aqju ,16162 + 2 ^ 0 ,j, 2 , 061 , (17) 


k l — k 
-1^2 


where (®2 ^^i/ 2 ('^)) 6.63) and bi^'^ is a Laplace coefficient (MD 6.67). We define Aj^n = 

O 2 A 0 J to simplify the expressions below. 

We rewrite TT.^) in complex notation (in the plane-parallel limit, expanded to second order in ei), giving: 


TZd =TZd,o + Re 


1 


( ^joo + Ajio6i + Ajoi 62 -I- 2^70262 + Ajiieie2 + -Aj2o6i) 


i>i 


1 


1 


'TZofl — 2 (^^000 + A01061 -I- A00162 + 2^00262 + Aoii6ie2 -I- -A02061J . 
where we have used the fact that Aoj^k,i = Ao-j^k,i since bi^'^ = bi~^\ Likewise, 

(l-kei)(l-f . 


(18) 


77 E = —Re 


(19) 
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Taking the derivative oiTZo and TZe with respect to 9i gives to first order in ei: 


— -^Sli — niiJ, 20 :Re f^joo + + ^ioie 2 ~ <a(l + ei ~ 2e2)(5ji) nje ^ 

tti UUi \ / 

Li>i 

where we have used nf = Grrii^/a^ = and fii = mijm-),. The derivative of TZi with respect to ei is: 


2 — n\n2aRe - ^Aqio + ^ 011^2 + ^02oei^ 

+ + 24^1162 + 24j2oei ~ Q:(1 — 2e2)Sji^ . 

j>i ^ 


For the outer planet, 


a| 862 ai ^ 


n^HiRe ^^ioo + ^jioei + ^joie2 ~ ot ^(1 + £2 — 2ei)6ji^ 

J>i ^ 


The derivative of 7^2 with respect to £2 is: 


^ —n^HiRe - ^diooi + ^oiiCi + ^002^2^ 


+ {^joi + ^iii£i + 24j02e2 — a ^(1 — 2ei)i5ji^ . 

j>i 


The angles and radii in the derivatives of the disturbing function can be expanded to first order in eccentricity, 
yielding: 


^= ni/r2a-Re V ((i^oo - a5ji) 

+ (j^ljoo ~ — ^adji)zie'^^ + ( ~ J^ioo ~ ^^jio + ^aSjijzie 

+ ( — jb 4 joo ~ + (j^ioo ~ ^^joi — 2aSji)z2e . 

c? ^ ni/X 2 Q:i?e 5^010 - 5 ^ 020216 “’^^ - 5 ^ 011 - 226 “’^'' + ^ | (i^io - aSji) 

1 L j>i \ 

+ {j^jw ~ 5^i20 ~ aSji^ jA^ig — ^Aj2o + Zic 

+ ~ ^Ajii^ + (jAjiQ — ^Ajii — 2a6ji^ Z2& • ( 24 ) 

where ^ = \\ — A2. We have also combined the j = 0 eccentricity terms since Re{z*e'^') = Re{zie~'^') = ^(z*e’^’ + 
For the outer planet, 

«2 092 ^ 

+ (j^ljoo ~ ^Ajio — 2a ^6ji)z^e^A + ( — jAjoo — 521^10)2:16 
+ ( — jAjoo — ^^Ijoi + fo! ^<5^1)2:26’^^ + (jbljoo ~ 5 ^ioi ~ 5Q! '^dji)z2e . 

\^^^=nlniRe 5^001 - 5^01121 - 5^002226*^" + ^ [ (Ajoi - a~'^5ji) 

«2 062 [ ^ V 

+ ^jbljoi — 5^jii ~ 2a ^Sji^ z*e^^^ + jbljoi — ^AjuJ zie 


jbljoi — 521^02 + Q; ^< 5 ^ 1 ^ 2:26*^^ 


^ 236*^^ + (^jAjoi - i 


~ 2^102 — Ct Sji ) 226 


-2 r . i ^-iAj 
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5.5. Trial solution 

The derivatives of the complex disturbing function contain terms that are proportional to for 

j > 0. We treat these terms as harmonic driving terms, and solve the inhomogeneous partial differential equations 
term by term. We expand the complex solutions for 66i and 6ei as trial solutions: 

j>i 

j>i 


fc=l,2 


d0ij=60^°]+ Y 


fc=l,2 


and 


6e2 = 6eiYeie'‘'^^~"''> + Y 

j>i 

S02 = S0[Yeie^^^^-^"'^ + Y 502, 

i>i 

5e2,j=5/°Y E '=)+<5e(7/^efee-‘(^'=—'=)) , 

k^l,2 

602,j=60^°] + Y . 


k^l,2 


We also define the solutions to zeroth order in eccentricity as: 


S0‘f'> =Yof^je'^'<’. 

j>i 

Then, the (real) transit timing variations are equal to 

Sti = —ni^Re{69i + 2eid9^'^) 


= 6t 


(- 2 ) 

1:0 

-1 i 


E 

i>i 


5t'ii 


E (^4: 


(+fc) 

j 


■6t[Y 


fc=l,2 


Stf]=-n^^Re{60^°]e^^*) 


= -n^ 


^Re [[50[Y - efee 




and 


5t2 = — ^-2 ^Re{502 + 2e2S02°^) 


= 6t 


(+ 1 ) 

2.0 


E 

i>i 


StfY-^2'M50^2,je'^^) 
5 tiY=-^ 2 "Re - < 54542 ) 


°^2,j 


E (->4^ 


(+fc) . ^A — k) 
3 ^^2,j 


fe=l,2 


(26) 


(27) 


(28) 


(29) 


(30) 


5.6. Inner planet coefficients 

Substituting the (j, 0, ±fc) trial solutions into the above differential equations, to zeroth order in eccentricity we find 
for the inner planet: 


2iffi 


4) 


(0)N 


—2if3j —{Pj+3)) 14^5/ V 4^10 — a4i 


ij(4joo — ctSji) 


(31) 


where Pj = j{ni — n 2 )/ni = j{l — a^/^) and we have divided the equations by nf to make dimensionless. 
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Similarly we can write down the equations for the coefficients to first order in ei. Note that to solve for we 

need to compute Hence we can subtract the matrix on the left times the vector from both 

sides of the equation. This may be rewritten in dimensionless form as: 


-{Pj ± 1)2 2i(/3, ± 1) 

-2i(/3,±l)-({^,±l)2 + 3), 




l±/3, -i(/3,±l)^ 
i(3/3j ±2) 5 , 


'^20: 


\j {±3 A 


yoo 2^110 


X 0 

+ i20i5ji{l =F 2)^'\ 


ijHjio 2^120 T OL8j\ 


(32) 


The solutions of equation 


31 


to solve for the first order eccentricity terms. 

To first order in 62 in dimensionless form (for j > 1), 


for (i50^°j, (5e^ j) may be plugged into the first term on the right hand side of this equation 


-n± 


2i?7± 

-2i77± -( 77 I + 3) 



= ^20! 




i ^Tjb 4 joo 2^101 (1 T l)Q!i 5 ji^ 

Tj^jio - ^Ajii - (It l)aSji ^ 


(33) 


where r]± = j3j ± — n- 2 )ln\ ± n 2 /ni = j(l — ± o?!'^ and Aju = —(2Ajio + Aj 2 o) = —‘2adh^^^^lda — 


'^d'^hfiJda^ 


— J\^ y_i_Li CA.XJ.V-1. I 

Note that for j = 1 the rj^ term becomes rj^ = 1. The frequency dependence of this term is at the 


Keplerian frequency of the inner planet, ni, and the determinant of the left hand matrix becomes zero as the second 
row of the matrix is equal to 2i times the first row. This singularity occurs due to the fact that the equations become 
those of a resonantly driven oscillator, which means that the amplitude grows linearly with time (or, equivalently on 
short timescales, the Keplerian frequency is shifted). This term is not relevant for transit-timing analyses as it occurs 
at the frequency of the transiting planet and grows on the secular timescale; consequently we will drop this term for 
now and discuss below in appendix B. 

For j = 0, 


/ _c3 _2ia3/2 ^ 




= M 2 Q! 


0 

-ii, 


on 


(34) 


5.7. Outer planet coefficients 

Defining Kj = jffii — n 2 )(n 2 = dividing this equation by gives the dimensionless form of: 


2iKi 


3 —••'3 

—2inj —(^2 T 3) 



and for the equations to first order in ei for j > 1 , 


-Cl 2ie± ' 

-2ie± -(Cl+ 3)2 


where C± = Kj ± a ^/ 2 ^ while for j = 0 , 


'seffi 


= Ml 


= Ml 


-ij{Ajoo - a "^Sji) 
Ajm — a~'^5ji 


-ij(±jHjoo - l^iio - (1 ± l)a ‘^5ji) 

^3Ajm - \Ajii - (1 ± l)a~'^Sji 


(35) 


(36) 


v-3 


2ia-3/2 


—2iia ^/2 



and for first order in 62 , 


— {Kj ± 1)2 2i(Kj ± 1) 

-2i(«,±l)-((«,±l)2T3) 


1 ± Kj —i(Kj ± 1) 

i(3«:j ±2) 5 ^ 



= Ml 


-p, 


on 


50(f) - 50(f 
(± 2 ) 


^ji^jAjOO 2 -^jOl ffi 2 (1^2)0 5^1) 

Tjf 01 - lAj02 ± a“2^^.^ 


(37) 


(38) 


The function u originates from the inversion of the matrices in the left hand side of equations (31)-(38l, which each 
have the form of: 

-72 2 ii 7 
^-2i7 -(72 T 3), 


(39) 






18 


Agol & Deck 


where 7 is the dimensionless frequency in units of the orbital frequency of the transiting planet. Since TTVs only 
depend upon 69, then only the first term in the inverse of this matrix times the right hand side of the equation gives 
the function u. These terms are driven by the disturbing function only. 

The function v± results from the driving ter ms caused by appearance of the zeroth-order eccentricity equation on 
the right hand side of the linearized equations (16). These can be expressed as the inverse of the zeroth-order matrix 
times the coefficients of the disturbing functionariving the TTVs at zeroth order; the first term on the right hand 
sides of equations ( [3^ and (36) times the inverse of the matrix on the left hand side yield the functions v±. 

The expressions for u and v± result from the coefficients in equations (M]|^) which can be found by inverting the 


matrices. Note that the coefficients of 69i are imaginary, while 5ei are real; thus, 59i, and hence Sti, will always have 
a sine dependence, while Set will always have a cosine dependence. 

As with the inner planet, for j = 1 the term becomes ^_ = 1. The frequency dependence of this term is at the 
Keplerian frequency of the outer planet, 712 , and hence the determinant of the left hand matrix becomes zero as the 
second row of the matrix is equal to 2ii times the first row (as occurs for the inner planet). We will drop this term for 
now and discuss next in appendix B. 

APPENDIX B 


SECULAR TERMS 

In the foregoing analysis we neglected the presence of secular terms, which in the disturbing function appear at zero 
frequency as well as at the Keplerian frequency of the planet that is being perturbed. These terms cause corrections 
of to the ephemeris of the planet, and thus cause an error of to the computation of a from the best-fit 

mean period. The correction to a affects the coefficients of the TTVs at order and so it can be neglected 

for the purposes of the first order in eccentricity transit timing solution. However, the solution we present here may 
have other applications, such as for radial-velocity planets or astrometric motion, which are not aliased at the orbital 
frequency of the planets, and so these secular terms enter at the 0{^^) level. In this appendix we compute these 
secular terms to hrst order in /i and e. 

In the equations of motion for the inner planet, to include the secular and Keplerian frequency terms, we will use 
angular momentum, li in lieu of angle 9i. The equations of motion become: 


h = 


dn^ 


ei - 




ai{l + e,r 


dOi ’ 

= -n?(l + ei )-2 


1 57^l 
af dei 


(40) 


where we have made the substitution ri = ai(l -|- ei) into equations (11) and we have divided by oi. In equation 24 
we keep only the secular terms and terms at frequency ni, giving: 


ei - 


n 


ii — —ril^ 2 aiaRe [ii(Aioo + 5^101)2:26’ 


,i Ai 


4 (l + ei)3 (1 + ei)^ 


-I- n\n 20 iRe 5A010 — 5^020^16’^^ — (Ahq -I- hA^ ''--*''*^1 


1111 


^ * i 


(41) 


for the inner planet, and similarly for the outer planet 


^2 — i(Aioo + 2 ^iio)- 2 :ie 


iAo 


^2 - 


Kl + e2)3 (l + e 2 )= 


+ \n2iJ-1Re Aooi — ^0022:26’'''^ — (2A101 -I- Aiii)z\e 




(42) 


(43) 


where we have taken the complex conjugate since this does not change the real component. 

The solution for the inner planet’s angular momentum to first order in eccentricity is: 

h = nia\{\ — |/i2aAoio — M2 Q:(Aioo + \AiQi)Re\z2e'^^^ . 

This can be substituted into the equation for ci, keeping terms of order eccentricity, to obtain: 

e'l -I- n\ei = -nlii2aRe[gie^^^] 

9i = 5(3Aoio + Ao2o)’2:i + ( 2 A 100 + Aioi + Ano -I- ^Aiii)z2- (44) 

Note that we have chosen the constant of integration in R to cancel the constant term in the disturbing function 
derivative that appears in the equation for ei so that the ei does not have an offset. This is because we prefer to 
specify the value of the semi-major axis in the initial conditions. 

The solution to this equation is: 


61 ,. 


= Re 


(- 2:1 -f i/i 2 Q;giAi/ 2 ) e 


i Ai 


(45) 


Note that this solution grows in amplitude linearly with time; however, the growth is slow, occuring on the secular 
timescale times the inverse of the eccentricity. 
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With these solutions in hand, we can solve for d\ = l\lr\ 
respect to time, giving: 


^(1 — 2ei) + 0(ef). We then integrate this with 


Oi^sec = M(l - \^X2aAow) + i?e|^ - 2izie’‘^^ - /i2Ma:gie’‘^^ 
— i / j , 20 : |(/i — (Aioo + 1^101)^2 ~ 5^010^11 e’'* 

A similar solution can be derived for the outer planet, with: 


I 2 = n2a^ — j/riAooi — /xi(Aioo + \AiiQ)Re ^ . 

As before, this can be substituted into the equation for € 2 '- 

€2 + nle 2 = -nlfXiRe[g 2 e^^^] 

32 = 5(3Aooi + ^002)^2 + (2A100 + Alio + Aioi + ^Aiii)zl. 


This equation has solution 


^ 2 ,sec 


Re 


{-Z 2 + 13231 ^ 2 / 2 ) 6 *^" 


(46) 

(47) 


(48) 

(49) 


Substituting this into the relation 62 = ~ 120 , 2 ^{^ — ‘^£ 2 )+ 0 ( 62 ) and integrating 62 with respect to time yields: 


^ 2 ,sec — ^ 2(1 — jMiAooi) + Re — 212 : 26 ’'''" — 3 iA 232 e’^" 

—131132 ~ (Aioo + 1 ^ 110)21 — ^A 02022 1 6*'''" 


(50) 


We have verified these solutions by plugging them back into the differential equations, both analytically and numerically. 
The solutions for 9i^sec can be transformed to timing variations with SU^sec = —9~g^^,{9i^sec — A^). 



